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We show that superconducting currents are generated around magnetic impurities and ferro¬ 
magnetic islands proximity coupled to superconductors with finite spin-orbit coupling. Using the 
Ginzburg-Landau theory, T-matrix calculation, as well as self-consistent numerical simulation on a 
lattice, we find a strong dependence of the current on the direction and magnitude of the magnetic 
moment. We establish that in the case of point magnetic impurities, the current is carried by the 
induced Yu-Shiba-Rusinov (YSR) subgap states. In the vicinity of the phase transition, where the 
YSR states cross at zero energy, the current increases dramatically. Furthermore, we show that the 
currents are orthogonal to the local spin polarization and, thus, can be probed by measuring the 
spin-polarized local density of states. 
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Superconductor-ferromagnet heterostructures were re¬ 
cently proposed as a viable platform for realizing topo¬ 
logical superconductivity (TS) [HE], which can host Ma¬ 
jor ana fermion quasiparticles at vortex cores and bound¬ 
aries mu. The Majorana fermions obey non-Abelian 
statistics and may be utilized for topological quantum 
computation [311]. The key ingredients driving these 
systems into the topologically nontrivial regime are the 
spin-orbit coupling (SOC) and magnetism. Recently, the 
search for experimental realizations of TS has also led to 
engineering the Yu-Shiba-Rusinov (YSR) [TQTIT2] bands 
induced by magnetic atoms on the surface of a supercon¬ 
ductor [T3f[25] . Following this recipe, zero-energy peaks 
in the tunneling spectrum were recently measured at the 
ends of a one-dimensional (ID) chain of magnetic atoms 
[26] . Such a tunneling spectrum could be the evidence of 
the Majorana edge states, although alternative explana¬ 
tions are also possible [27] . 

The interplay of SOC and magnetism has another re¬ 
markable consequence. Consider a two-dimensional (2D) 

surface of a three-dimensional (3D) material. The effec- 

2 

tive Hamiltonian of the surface h(p) = + A (cr x p) z 

contains the Rashba SOC due to the absence of the inver¬ 
sion symmetry at the surface. Then, the velocity opera¬ 
tor v = = ^ -\- X z x cr contains a spin-dependent 

term that gives an extra contribution to the current 

Jextra ~ Xz X (<t) . (1) 

A ferromagnet proximity coupled to the superconductor 
would render a finite spin polarization (cr) ^ 0 and, thus, 
generate a current as schematically shown in Fig. [lja). 
The phenomenon of driving a current with magnetism 
is known as the magnetoelectric effect. This effect may 
vanish in metals due to dissipation but survives in su¬ 
perconductors lacking inversion symmetry [281132] . The 
magnetoelectric effect was also recently discussed in a 
pure ID model of TS [ 33] . 
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FIG. 1. (color online) (a) Schematic representation of the non¬ 
local currents (red arrows) induced by a ferromagnetic (FM) 
island on the surface of a superconductor with the Rashba 
SOC. GL solutions for the current around a circular ferro¬ 
magnetic island (gray area) with S — S z (b) and S = S x 

(c). 


In this work, we show that the magnetoelectric current 
is universally generated around single magnetic impuri¬ 
ties and ferromagnetic islands, which have been recently 
studied in the context of TS [T31I26] . More specifically, we 
first derive the extra terms in the Ginzburg-Landau (GL) 
free energy corresponding to Eq. 0- For a small ferro¬ 
magnetic island on a superconductor with SOC, we find a 
strong dependence of the current on the relative orienta¬ 
tion of the ferromagnetic moment. The current circulates 
around the ferromagnetic island and is short ranged for 
the ferromagnetic moment normal to the surface. On 
the other hand, the current has a dipolar power law de¬ 
cay for the ferromagnetic moment parallel to the surface. 
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Next, we discuss the current generated around a point 
magnetic impurity and show that the current is carried 
by the impurity-induced YSR states. We also perform 
a self-consistent numerical calculation and find a strong 
nonmonotonic dependence of the current on the strength 
of the ferromagnetic moment. The current strongly peaks 
at the phase transition, where the YSR states cross zero 
energy E = 0. We further demonstrate that the cur¬ 
rent can be mapped by measuring the spin-polarized lo¬ 
cal density of states (SP LDOS), which, thus, provides 
a probe of both the current and the phase transition. 
Our findings are, therefore, highly relevant for the ongo¬ 
ing search of the Majorana bound states in ferromagnetic 
chains [T3ti26] . 

Ginzburg-Landau treatment. - We start by considering 
a ferromagnetic island deposited on a 2D surface of a con¬ 
ventional s-wave superconductor with the Rashba SOC 
as illustrated in Fig. [lja) and described by the Hamilto¬ 
nian 

H = I J d 2 r \l/ t (r) [h(p)r z + A t x - S(r) • er] ’F(r), 

h{p) = — + A (<r x p) z - p, p = —i(V x , V y ). (2) 

Here ’J = (ip- f, ipi, +j, —+j) T is a four component spinor, 
a and r are the Pauli matrices acting in the spin and 
particle-hole Nambu space, A is the superconducting 
gap, and we set e = h = 1. The ferromagnet and its 
coupling to the superconductor are described by the spa¬ 
tially dependent vector S = (S x , S y , S z ). An intuitive 
and qualitatively correct picture of the currents can be 
derived using the GL free energy 


F = d 2 r 


+a(z x S)-A+P(VS Z x z)-A 

V6> 


A — (Ax, Ay) — A H——, 


(3) 


which is valid at length scales larger than the supercon¬ 
ducting coherence length £ sc . In the first term propor¬ 
tional to the superfluid density n s , vector A encapsu¬ 
lates both the superconducting phase 0 and the vector 
potential A. The second and third terms describe the 
coupling between the Rashba SOC and magnetism and 
are derived in the appendix. For example, in the limit 
Pf A m\ 2 A, the coefficients are [34] a = mA/27r, 
(3 = m 2 A 2 /, and, thus, only present at finite SOC, 
i.e. when A / 0. The term proportional to <a, known as 
the magnetoelectric term ( 2 HH 201 E 21 HU ESI, is allowed 
only in the absence of inversion symmetry. 

Within the above framework, we now discuss the cur¬ 
rents induced by a ferromagnetic island of a uniform 
disc geometry, which we model as S(r) = S Oh (R — r), 
where Qh(z) is the Heaviside theta function. We find the 
current from Eq. § as 
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SF 

Ja 


A =0 


n 

-^V0 + a(zxS)+(3(VS z xz). (4) 


First consider the an out-of-plane ferromagnetic moment 
S = Sz and 0 constant. Then the current is given 
by the last term in Eq. ©• The current is localized 
near the boundary as j(r ) = —/3S(r x z)S(r — R ) and 
circulates around the ferromagnetic island as shown in 
Fig. gb). Since the GL equations are valid at r > £ sc , 
the 6 function in the current solution is artificially broad¬ 
ened to a scale of the superconducting coherence length 
£ sc for visualization purposes. For an in-plane moment 
S = Sx, both of the first two terms in Eq. © are 
nonzero. The contribution given by the term a(z x S ) 
is constant over the region covered by the ferromagnetic 
region and discontinuous at the boundary. However, the 
first term ^ V# fixes this discontinuity. Indeed, the vari¬ 
ation of the free energy over 0 gives the continuity equa¬ 
tion: 0 = V • j = + a V • (z x S). The last 

expression is the 2D Poisson equation with a source term 
that we solve for 0 and plot the currents in Fig.[ljc), see 
appendix for more details. The current is constant over 
the region covered by the ferromagnet, i-e. j(r) = jgf 
for r < R, and has a dipolar profile outside of it, i.e. 
j (r) = 2r ( 4 "^ — ^ for r > R. Here the effective “dipole” 
moment is defined d = aR 2 ^(zxS). We also note that if 
the coefficients a and /3 are large, vortex solutions for the 
superconducting phase 0 are favored by the free energy 
expression Eq. ©• 

Microscopic calculation.- To complement the above 
GL analysis, we also study microscopically the cur¬ 
rents generated around a single point magnetic impurity, 
i.e. we set S = S S(r) in the Hamiltonian Eq. ©• In con¬ 
trast to the GL approach, the Green’s function method, 
used below, allows to study effects to infinite order in 
S and also at distances smaller than the superconduct¬ 
ing coherence length, i.e. for r <C £ sc - We evaluate the 
Green’s function of the superconductor in the T-matrix 
approximation 


G rr '(u) =3W'M +5 r oMT(w)£<>'M, 




-S (T 

1 + 5 • ergooM' 


(5) 

(6) 


The Green’s function of a clean superconductor in real 
space at r is (for r <C £sc) 


g r0 (uj) = -ir-^=== [fo(r) + i(<r x r) /i(r)], 


(7) 


where f 0 (r) = \ [p + Jo(Ppr) + p J 0 (p F r )] and fi(r) = 
^ [p + Ji{Pf t ) ~ P~ Jo and J\ are Bessel func¬ 

tions, p p ~ pf T Am are the Fermi momenta of the spin- 

polarized Rashba bands, and p ± « are the corre¬ 

sponding density of states (po = m/ir). Equation 0 is 
calculated with the assumption p A > 0. The second, 
spin-dependent, term in Eq. 0 is a consequence of the 
Rashba SOC and vanishes if A = 0 . The poles of the 
T-matrix give the energies of the impurity-induced YSR 
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FIG. 2. (color online) (a),(c): Currents around a point magnetic impurity calculated using the T-matrix approach. The 
impurity is located at ro = 0, and the direction of the magnetic moment S is shown in the inset. The vectors in green indicate 
the direction of the in-plane spin polarization determined from SP LDOS. The current and the spin polarization are orthogonal, 
which is consistent with Eq. |l]). In order to enhance the figures, the currents are plotted away from the impurity, i.e. for 
Pft > 2.5. (b),(d): SP LDOS calculated at n. The insets show the SP-LDOS in the vicinity of the positive YSR state. 
Red and blue lines correspond to SP-LDOS calculated for opposite directions x and —x, whereas the dashed green line is the 
resulting spin polarization. Panels (a) and (b) correspond to the out-of-plane magnetic moment, whereas panels (c) and (d) 
correspond to the in-plane magnetic moment. 


subgap states IMHI M 
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which are unaffected by the Rashba SOC [35]. The ener¬ 
gies of the YSR states, however, depend on the ferromag¬ 
netic vector magnitude S. For a critical value S = 2/7Tp, 
the energies of the YSR states reach E = 0, and the 
system undergoes a quantum phase transition as the two 
YSR states cross [Hamm. For simplicity, let us temporary 
fix S' = 2/V^Trp, which corresponds to Ey SR = d=A/2. 

The current is equal to the expectation value of the ve¬ 
locity operator, which can be expressed using the Green’s 
function as 


o( r ) 


lim 

r' —>r 
<5—^+0 


/ 


duo e iujS 
2lH 


x Tr 


1 +T Z 
2 



+ A z x cr 


(9) 


G < pjt 1 ( uj ) 


r ^ 1 /pf- Note that the current in panel (c) is not con¬ 
tinuous. This can be understood by using the analogy 
with the Ginzburg-Landau current Q. For the in-plane 
vector S', the current consists of the bare term a(z x S), 
as well as the condensate term ^ V0. These two distinct 
contributions to the current are discontinuous, however, 
their sum is continuous. Since, the T-matrix calculation 
is not self-consistent, it does not take into account the 
reaction of the condensate that would fix the discontinu¬ 
ity. We discuss a fully self-consistent calculation, which 
demonstrates the continuity of the currents, in the next 
section, as well as in appendix. 

According to Eq. 0 the current and the spin polar¬ 
ization are coupled. Thus, we expect a nonzero in-plane 
spin polarization even away from the impurity site that 
sustains the nonlocal currents shown in Figs. |2ja) and 
(c). So, we evaluate the SP LDOS using the Green’s 
function as 




— —ImTr 

IT 


1 ~\~ T z 1 + G j 
~2 2 


G rr (uj + id) , 


( 10 ) 


In addition to the usual gradient term [401 in the paren¬ 
thesis, there is also a spin-dependent contribution due to 
the Rashba SOC. We evaluate the current in Eq. Q using 
the Green’s function in Eq. 0 and plot it in Figs. 0a) 
and (c) for the cases of out-of-plane S = Sz and in¬ 
plane S = Sx moments, respectively. We note that only 
the pole in the T-matrix corresponding to the YSR state 
gives rise to nonzero currents. Both panels show con¬ 
centric patterns of current centered around the impurity. 
In the case of the out-of-plane moment, the current cir¬ 
culates around the impurity. In contrast, in the case 
of the in-plane moment orientation, the current points 
predominantly in the y direction. The currents shown 
in Fig. [2|a) and (c) for point magnetic impurities agree 
qualitatively with the patterns obtained for the circu¬ 
lar island within the GL theory and shown in Figs, [ljb) 
and (c). However, in contrast with Fig. [l] the currents 
in Fig. [2] display fine Friedel oscillations on the scale of 


where j = x,y,z denotes the polarization axis. From the 
SP LDOS we also define the energy-dependent local spin 
polarization as 

= -/v'M- (n) 

In Figs. [0 (b) and (d), we plot both the SP LDOS and the 
spin polarization at the point rq with solid and dashed 
lines, respectively. First, consider the out-of-plane mag¬ 
netic moment S = Sz in panel (b). The SP LDOS 
peaks at the superconducting coherence peak, i.e. at 
uj = A, as well as at the subgap YSR state energy, 
i.e. at uj = £Asr = A/2 . The SP LDOS correspond¬ 
ing to the opposite directions j = ±x, shown with red 
and blue lines, are notably different at the YSR state. 
Therefore, the YSR state has a finite spin polarization 
along the x axis, shown with a dashed green line. This 
feature is a consequence of the spin structure of the 
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FIG. 3. (color online) Self-consistent numerical calculation of the currents induced by a point magnetic impurity with the 
out-of-plane moment S = S z. (a) Bogolyubov-de Gennes spectrum (top) and magnitude of current j (bottom) as a function 
of S. (b)-(d) Spatial profile of the currents plotted at discrete points on the lattice for increasing magnitude of S as indicated 
in panel (a). The current reaches the maximum value j c and switches direction for S = S c , i.e. where the YSR states cross 
zero energy. Note that in order to enhance visibility of the current in panels (b) and (d), the arrows representing the current 
are magnified tenfold as indicated by the magnification ratio in the top-left corner of each panel. 


Green’s function Eq. 0 and vanishes in the absence of 
the Rashba SOC. Now, consider panel (d) correspond¬ 
ing to an in-plane moment S = Sx. The YSR state 
in this case has a dominating spin polarization in the 
+x direction with only a small admixture of the oppo¬ 
site spin. In Figs. |2^a) and (c), we plot the direction of 
the in-plane spin polarization for the positive YSR state 
a r (Ey SR ) = [<(£+ sr ) ct^(Ey SR )] at the point r 0 = 0 
as well as rq = 7x/pp and r 2 = 7y/pF- Note that the 
spin polarization of the negative YSR state is opposite, 
i.e. cr r (Ey SR ) = —a r (Ey SR ). The current and spin po¬ 
larization are consistently orthogonal, which agrees with 
Eq. Q. So, it is possible to map the current generated 
by magnetic impurities and ferromagnetic islands using 
spin-polarized scanning tunneling microscopy (SP STM). 

Self-consistent numerical modeling.- The T-matrix 
approximation discussed above predicts currents which 
are qualitatively consistent with the GL results. How¬ 
ever, the T-matrix approach does not capture the influ¬ 
ence of the magnetic impurity on the superconducting 
order parameter. It is known that the superconducting 
order is strongly renormalized and may even change sign 
[371131 SB in the vicinity of the magnetic impurity. In or¬ 
der to take this into account, we also perform a fully self- 
consistent numerical simulation [[42] of the point magnetic 
impurity on a lattice [43ll45l and show the results for an 
out-of-plane magnetic moment S = Sz in Fig. [3] Panels 
(b)-(d) show the current for increasing values of the fer¬ 
romagnetic moment S. Note that the Friedel oscillations 
are not fully visible here since the calculation is done for 
a coherence length such that £ sc < 1/pf- In panel (a) 
we show the Bogolyubov-de Gennes spectrum (top) and 
the magnitude of the current (bottom) as a function of 
S. For small S (b), the current circles around the impu¬ 
rity, which is consistent with both previous Figs, [l] and 
[2] With further increase of S', the current grows and ul¬ 
timately undergoes a first-order discontinuous transition 
at a critical value of magnetic vector S = S c . There, the 
current abruptly reverses direction, as shown in Fig.J^c) 
and (d), and reaches its maximal magnitude. This is 
accompanied by the YSR states crossing at zero energy, 


and the superconducting order parameter reversing sign 
at the impurity site. We note that the YSR states also 
have a first-order avoided crossing at zero energy [39] as 
shown in Fig.[3ja). With further increase of S, supercon¬ 
ductivity is suppressed and the currents diminish in the 
vicinity of the impurity. More details on the numerical 
simulation can be found in the appendix. 

Concluding remarks.- We have shown that super¬ 
conducting currents are generated by ferromagnetic is¬ 
lands and single magnetic impurities in 2D supercon¬ 
ductors with spin-orbit coupling. The currents originate 
from the magnetoelectric effect and are a direct conse¬ 
quence of combining SOC and magnetism. The discussed 
currents are unavoidable in ferromagnet-superconductor 
heterostructures, which have been proposed as a plat¬ 
form for topological superconductivity with the Majo- 
rana boundary states HH3III3H23I. We find a strong de¬ 
pendence of both the spatial pattern and magnitude of 
the currents on the direction of the ferromagnetic mo¬ 
ment. The currents are localized on the scale of the co¬ 
herence length in the case of the out-of-plane local mag¬ 
netic moment, whereas the currents have a dipolar power 
law decay in the case of the in-plane magnetic moment. 
The presence of these non-local currents may induce long- 
range interactions between local magnetic moments on a 
superconductor [46], which could qualitatively change the 
behavior of the Majorana modes in such systems. Fur¬ 
thermore, by analyzing the currents in detail, we find 
that they are carried by the subgap YSR states induced 
by point magnetic impurities. The YSR states are spin- 
polarized, and the current is orthogonal to the local spin 
polarization. Moreover, the current magnitude peaks 
sharply at the phase transition, where the YSR states 
cross at zero energy. Thus, by using SP STM it should 
be possible to map out the currents as well as detect the 
phase transition, which is paramount for finding TS and 
the Majorana modes. 
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Appendix A: Derivation of the extra terms in the 
Ginzburg-Landau theory 


The non-local coupling of the vector *4, = (A x ,A y ) = 
A + and ferromagnetic polarization S in 2D momen¬ 
tum q = (q x , q y ) space is given by the extra term in the 
free energy 


The tensor K ab 
fermions 


is obtained by integrating out the 


K ab {q) = f | Tr 

UJ,p 


Va,g p+ 1 (iu) v b g p _ a (• iu ) 


(A2) 


where the integration variables are given below the inte¬ 
gral sign for brevity, i.e. f = f The integration 

u.p 

over continuous frequency uj corresponds to zero temper¬ 
ature T = 0. The factor 1/2 takes care of the doubling of 
degrees of freedom in the 4-by-4 Bogolyubov-de Gennes 
representation. The Green’s function as well as the ve¬ 
locity operators are defined as 

9P{IU) = ^-[£CP)-A(<TX P )jT,-ATy (A3) 

v = — + A (z x cr). (A4) 

m 

We expand the tensor K ab (q) up to the first order in q 

K ab (q ) = K ah { 0) + q c d q K ab { 0) + 0(q 2 ). (A5) 

and find the tensors K ab ( 0) and d c K ab ( 0) in the following 
two subsections 


T-afr(O) ^ zba ^ 0 5 (A6) 

Z7T 

2 a 2 

d c K ab (0)=i6 az e zcb p, P= T ^. (A7) 

47T Pp 

Here, the coefficients a and /3 were evaluated in the limit 
of the small superconducting gap A and large Fermi mo¬ 
mentum pf = ^m/x, i.e. A pp rn\ 2 A. It may 
also be useful to instead express the coefficients using 
the density of states po = W 71 " an d Rashba momentum 
Pr = mX as 


« = P 0 a1, /? = A,a||. (A8) 


We can now sub stitu te Eqs. (A6) and ( |A7| ) in Eq. (A5) 
and rewrite Eq. (Al) in real space r = (x~y) as 

Text r a — J a(zxS)-A + /3(VS z xz)-A. (A9) 


The term a was also recently derived in Ref. [36]. We 
are not aware of a previous derivation of the term /3. 


1. Derivation of the tensor K ab ( 0). 

Let us calculate the first term in this expansion, i.e. 


Text r a — / S a (q)K ab (q)A b (-q). (Al) 

Q 


Kab{ o) = \ J Tr Wa, g P ( P ) v b g p (p)} . (A10) 

UJ,p 
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We define the operators 

n± = \ [i± (o- x p) z ], p = (Aii) 

Z p 

which project onto the spin eigenstates of the Rashba 
coupling A (<t x p) z corresponding to the two eigenstates 
=bA p. We expa nd the identity operator I via the projec¬ 
tion operators ( |A 11 | ) as I = II + + II_ and simplify the 
Green’s function 

= 14 {iuf-EUp) +11 - (,^-g(p) ' 

(A 12 ) 

where we defined the energies due to the Rashba splitting 


f±(p) = C(p) ± M #±(p) = \/C±(p) + A 2 . (A13) 


We substitute Eq. (A 12 ) in Eq. ( |A 10 | ) and obtain four 
terms. Only the terms containing both II + and II_ have 
distinct poles, which produce a non-vanishing result upon 
the frequency integration. Thus we obtain 

K ab (0) = ^ J Tr<r {a a U + v b n.+a a U_v b U + } 


/ 


Tr T {[iw + C+(p)r 3 + Ari][iw + £_(p)r 3 + An]} 

[{iuy-Ei{p)}[{iuy-Epp)} 


(A14) 

where the first trace is over the spin space, whereas sec¬ 
ond trace is over the Nambu space. The trace in the sec¬ 
ond line of Eq. (A14) gives 2[— uo 2 + £ + (p)£ - (p) + A 2 ]. 
Using the angular integration in momentum space, we 
simplify the trace over spin matrices in the first line of 
Eq. (A14) to Xe za i). Then, we integrate over frequency 


and obtain 


AT a fr( 0 ) — tzba 2 


A f E + (p)E-(p) - S + (p)£-(p) - A 2 


E + (p)E-(p)[E + (p) + E-(p)] * 

(A15) 


We can evaluate the integral in various limiting cases 

zv rn\ _ m ^ f 2A 2 p|,/3A 2 , p > A » Api? > raA 2 , 
ab(Uj “ 6zba 2 tt 1 1, Ap F > mA 2 > A, 

(A16) 


where pp = y/2m p is the Fermi momentum. The first 
line of Eq. (A16) can also be obtained from Eq. (7) of 
Ref. l36l. 


2. Derivation of the tensor d qc K a b( 0). 


We expand Eq. (A 2 ) and obtain the first order coeffi¬ 
cient 

dq c K ab = i J TV {< 7 0 [d Pc g p (iuj) v b g p (iu) 

oo,p 

vb ^Pcdpi^^W • 


We integrate by parts to shift the position of the derivate 
9 Pc , use t he ide ntity d Pc g = gr z v c g (which follows from 


Eqs. (A3) and ( |A4| )) and obtain 
&q c Kab = 2 j dpcdptyw) Vi) 


00, p 


-§ / TV {cr a gp(iw)r z i; c gp(iw) v 6 5 P (iw)} 

00,p 

Below, we omit additional terms containing powers of 
the Green’s function g n under the trace because they 
vanish upon the frequency integration (since poles lie on 
the same side of the imaginary plane). We substitute the 
expressions for the velocity Eq. (A4), omit the terms that 


vanish under the angular integration in the momentum 
space and simplify the equation to 


dq c K ab ( 0) = 


^ 2 tzcct zbb 


J Tr {a a g p (iu) r z a d g p (iuj) a b g p (iuj)} 


(A17) 

Here, we dropped the term containing 
PijP c Ti{a a g p (iuj)r z gp(iuj)} which vanishes upon the 
frequency integration. W e sub stitute the ex pansion of 
the Green’s functions Eq. ( |A 12 ) in Eq. (A17) and traces 
over the Nambu and spin matrices decouple such that 


E 


Ttv {<t 0 n. Sl < 7 g n S2 a,- n S3 } 


Sl ,, 2 , S 3=±[( iw ) 2 - £ S 2 »HM 2 - £ s 2 2 (p)HM 2 - eKp)} 

xTr T {[ico + £ Si (p)t z + At x ]t z [iuj + £ S2 (p)t z + A t x ] 
x [ico + £ S3 (p)r z + At x }} . 

First we evaluate the trace over the spin matrices 

t v CT {cr a n Sl a 5 n S 2 cr s n S3 } -y (1 - sis 3 ) 


where the angular integration in the momentum space 
was invoked to simplify the expression. The remaining 
trace over the Nambu matrices can also be evaluated as 

Tty {[iui + £ Sl (p)r z + AT x ]r z [iu + £ S2 (p)t z + At x ] 

x N + £,s 3 (p)t z + At*]} = -2J 2 K Sl (p) + 6 2 (p) + 63 (p)] 
+ 26i (p) 6 2 (p)6 3 (p) + 2 A 2 [61 (p) + 62 (p) - 63 (p)] • 

We integ rate over frequency, substitute all terms in 
Eq. ( |A17 ) and obtain 

dq c Kab{ 0 ) = (A18) 


^az^zcb f XA 


/ 


£ 2 (p)p 


C(p)Ap C(p)Ap _L_ 1 

^(p) ^(p) E+(p) E.{p) 


Similar to Eq. (A16), we evaluate the integral in the limit 
of small A and obtain 


dq c Kabify — ^az^zcb 


2 A 2 


47TP' 


2 ’ 


Xpf mX 2 A. (A19) 







































Appendix B: Ginzburg-Landau solution of currents 
induced by a ferromagnetic disc 

In this section, we provide details of the calculation of 
the current 


Ja 


A =0 


= 7 r L V0 + a(zxS) + ]3(zxVS z ). (Bl) 

2m 


around a ferromagnetic region with disc geometry 

S(r) = Se H (R-r), (B2) 

where the index in Oh denotes the Heaviside theta func¬ 
tion to contrast it with the phase of the condensate 0. 
Let us first consider the case where the spin is out-of¬ 


plan e, i.e . S = Sz. Then, taking into account Eq. (B2), 
Eq. ( |B1| ) becomes 


j(r) = ^V0^p[zxVS z {r)} 
= /3S (f x z) S(r — R). 


o=o 


(B3) 

(B4) 


The current is thus localized around the boundary of the 
ferromagnetic region as shown in Fig. [ljb) . Note that 
Eq. (B4) corresponds to vanishing superconducting phase 


0 = 0, which is valid for small constant /3. For larger val¬ 
ues of /3, the vortex configuration of the superconducting 
phase, i.e. 0(r) = atan (y/x ), m inimizes the free energy 
and the full expression Eq. ( |B3| ) for the current must be 
used. 

Now let us consider the case of an in-plane ferromag¬ 
netic vector, i.e. S = Sx. Then, the current Eq. (Bl) 
becomes 




(B5) 


Notice here that the Euler-Lagrange equation for the su¬ 
perconducting phase 6 gives the continuity equation for 
the current 


n s 


0 = V • j = AO + a V • [z x S(r )]. 
2 m 


(B6) 


(a) 


(b) 


A r ~ l/p F 



j l -y 


X 


X 


Ar ~ 1/pf 


t / 

t t 
t t 
\ \ 


i T i 


/ 

\ / 


t / A+ 

t f / 

\ \ t 

\ \ 

t t 

/ / 

/ t 


\ 


t t 

Sy ' 

T X 


t t \ 

t \ \ 

I N \ 

l * \ 
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FIG. 4. (Color online.) Current around a point magnetic 
impurity in a square tight-binding model obtained using self- 
consistent calculations. The lattice constant is a and the cur¬ 
rent vector is indicated on each tight-binding site with the 
impurity marked with a black dot. The direction of the im¬ 
purity magnetization is shown in the insets. 


We rewrite Eq. (B6) as AO = z • [V x S(r)] = 
— z • (f x S)S(r — R ), use the Green’s function 
Gl(t) = 2 ^: ln(r) for the 2D Laplace operator, which 
satisfies A Gl(t) = £ 2 (r*), and find the solution 
0(r) = - 2mac » t5oi f d 2 r' G L (r - r')r' a 5{r' - R), which 
after integration gives a simple result 


where the “dipole” 
d = aR 2 (z x S)/2 . 


moment is defined 


as 


e(r) = 


ma 

2 n s r 2 

ma 
- r 

n» 


r • (S x z) [R 2 +r 2 -\R 2 -' 


(Sxz) 


J 1 ’ 

i & 

L r 2 ’ 


r < R, 
r>R. 


Using Eqs. (BT) and (B2) we calculate the current (B5) 


j(r) = 


d 

R 2 J , 

2 r ( dr ) 


r < R, 
- 4 , r > R, 


(B7) Note that the problem of calculating the current dis¬ 

cussed above is formally equivalent to the magnetostatics 
(B8) problem of calculating the magnetic field B — H + 47 tAT 
induced by the ferromagnet of magnetization M. In 
the magnetostatics problem, the constant magnetization 
M is equivalent to a [z x S(r)], and t he m agnetic field 
strength H is the term in Eq. (B5). The diver- 

(B9) genceless magnetic induction B is equivalent to the cur¬ 
rent j. 
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Appendix C: T-matrix calculation 


which after insertion of projectors (1 ±r x )/2 becomes 


1. Green’s function in the real space 


We calculate the Green’s function in the real space as 


g r0 (u) = j 


i pr 


U - K(p) + -M<t x p)z]t z - A T x 


(Cl) 


We substitute the expansion of the Green’s function 
(fAl2|) in Eq. ((Cl]) 


9r o(w) = ^ J 


+ 


1 + i(cr x p) 
u - C+(p)t* - A t x 

1 — i(cr x p) 


~ A r x 


where £±(p) = £(p) =t Ap. The angular integration in 
the momentum space transforms the Rashba term in the 
momentum space (cr x p) into a corresponding term in 
real space and produces the Bessel functions Jo and J\ 


9ro(u) = ^J 


+ 


Jo(pr) + i(cr x r) Ji(pr) 

W - £+0) r * - Ar X 
Jo(pr) — i(cr x r) J\(jpr) 


u ~ f-( p)r z - A r x 


For small distance r C we can substitute the momenta 
in the Bessel function with their average values p = pp, 
which satisfy equations £± (p^) = 0. Then numerators 
do not depend on the integration variable p, and, thus, 
integration gives 

flroM = - 71- I/o( r ) + i ( (T X r ) /l( r )l - ( C2 ) 


v / A^' 


ClT 


where f 0 (r) = - [p+J 0 (pp) + p-J 0 (p F r)\ , 

fi(r) = 7^ [p+MpP) ~ P-MPf t )] ■ 


In Eq. (C2), the density of states for each branch of the 


Rashba spectrum has been defined as 


P± 


Po 


± ’ 


2 l±m\/pp 


Po = 


m 


(C3) 


For example at r = 0, the Green’s function has a conven¬ 
tional form 


9 oo M = - 
where p = p + + p_. 


irp uo + A r x 


2 VA 2 - cj 2 ’ 


(C4) 


2. T-matrix and spin polarization 

Using Eq. ( |C4| ), we can expand the T-matrix in Eq. © 

-Scr 


2 » = 


1 - 


TT£ S (T u-l-Ar, ’ 
2 VA 2 -w 2 


1 + STj; 


-Scr 


S— i 


Trp S cr uj+sA 
2 VA 2 -w 2 


We multiply the numerator and denominator of the frac¬ 
tion by 1 + 2 e p m i na te the cr matrices in 

the denominator 


E l + 8T X 
2 


—5 • < 


( 7 rpS) 2 (oj+sA) 2 
4 A 2 -w 2 


irp S ■ cr 


■ sA 


V^A 2 


■ CJ" 


and after a few transformations rewrite the expression in 
the final form 


= £ 


1 + ST X 


sA — i 


S=± 


1 + 


(irpS) 2 


(CJ JAysRy 


S-tr¬ 


ips 2 


L0 - 


E* 
■ sA 


VA 2 - cj 2 


(C5) 


where the YSR energies Ey SR are given in Eq. (pi). We 


substitute Eq. (C5) in the expression for the current © 
and find that the last term in the parenthesis on the sec¬ 
ond line vanishes after taking the trace. The remaining 
term proportional to S • cr has a pole at the YSR energy 
and gives a non-zero contribution to the current. Note 
that, in general, we expect the contribution to the current 
both from the localized subgap states and the delocalized 
supragap states. 


Using the expansion of the T-matrix in Eq. (C5), we 
calculate the spin-polarized LDOS (10) in the vicinity of 
the positive YSR state 


Pr*>) = [/o(0±r/i(r)] 


2 SA5(u> Ey SR ) 


1 + 


(~2~ )' 


-1 2 


(C6) 


for, e.g., perpendicular local moment S = Sz. The first 
term in square brackets is responsible for distinct SP- 
LDOS in the oppos ite directions ±®, whereas the sec¬ 
ond fraction in Eq. (C6) determines the overall strength 
of the YSR state. The terms /o(r) and r/i(r), defined 
in Eq. (C2), are of the same order sufficiently far from 


the impurity, and the YSR state acquires strong in-plane 
spin polarization = p+? - p~? « p+? according to 
Eq. (C6). Such a large spin polarization should be pos¬ 


sible to detect experimentally using SP-STM. 
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Appendix D: Numerical simulation 


tor Pi = J2ia C L C ™ ( see e -g- Ref - S2) : 


We consider the following tight-binding model on a 
square lattice (see e.g. Ref. 005]): 


dpi 

dt 



Let a and a be opposite spins and define 


(D3) 


H — -f^kin + H so + H sc (Dl) 

^kin = c\ a c jcr ~ M C icr C icn 

(i,j),cr i,c t 

Hs0 = y C UbaZ • X b W Cia' + H.C., 

ib 

Hsc = ^ + He.. 

i 

-f^imp = ~ c Q a (5 • 0’) aa , CQ a ' , 


where i and j are site indices and b is a unit vector point¬ 
ing along one of the four types of bonds on the square 
lattice. The parameters of the model are the strength of 
the nearest neighbor hopping t, the chemical potential /i, 
and the ferromagnetic vector S. Also, a superconducting 
pair potential v sc is used to self-consistently determine 
the superconducting order parameter through 


— ^sc 

= ~ V SC V ^i Ui t- ( D2 ) 

E v < 0 


We solve Eqs. and ( |D2[ ) self-consistently and cal¬ 

culate currents around point magnetic impurities using 
the expressions derived below in Appendix |D 1| Figure [3] 
corresponds to the following parameters: t = 1, fi = —4, 
A = 0.56, and v sc = 5.36. These values give a generic 
band structure of a lightly hole-doped Rashba SOC semi¬ 
conductor in proximity to a conventional 8-wave super¬ 
conductor. Panels (b), (c) and (d) correspond to a mag¬ 
netic impurity with S = 1.6,2.72, and 7.36 respectively. 
Figure 4 is plotted for the parameters t = 1, /x = —4, 
A = 0.56, v sc = 5.36, and S = 2.56, which is just below 

5 C . 


~~ Ci aa C ia C ’ 10 ' ' a aa C ia Cicr ’ 


J(j(j — ^ a crcr C i+bcr Cicr + a <7Cr C i a Ci + bcr ) 5 

J% = - (4 b a4+b«x C i<7 + 4 b a 4* C i+ba) , (D4) 


where b runs over the vectors from site i to its nearest 
neighbors and a ^ a , a ^, and a ^ are coefficients to be de¬ 


termined. It can be shown that Eq. (D3) can be written 
as 


cr b a 


(D5) 


It is clear from the definition of Eq. (D4) that S^ is an 


on-site source operator which converts a spins into cr, 
while is a current operator for cr-spins away from 
site i along bond b. Similarly is a current operator 
for cr-spins away from site i along bond b, but for spins 
which are flipped in the transfer process. Summing the 
two types of currents for both spin species, we arrive at 
the expression 


J ib =^j£„ (D6) 

oo' 

for the total operator for currents away from site i along 
bond b. The current operator on site i can now be defined 
as 


J i =^b J ib . (D7) 

b 


Evaluation of the commutator in Eq. (D3) reveals that 
the relevant coefficients in Eq. (D4) are 


4 b = - 

4 b = - £ • (o- x b)- a 1 (D8) 


Using Eq. (D4), (D7), and (D8) the current can finally 
be calculated as 


1. Calculating current 

Numerical expressions for the current are derived by 
considering the time rate of change of the density opera- 


^ b {it ('^i+ba'^io’ ^ia^i+bcr) 

h<j,E u <0 

+ 4 [z- (ax b)- a Vi +hff u ia 

+ z-(ax b)*- a vf a u i+bd -\} . (D9) 










